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Abstract 

We analyze the experimental hadron yield ratios for central nucleus-nucleus colli- 
sions in terms of thermal model calculations over a broad energy range, ^/snJj=2.7- 
200 GeV. The fits of the experimental data with the model calculations provide the 
thermal parameters, temperature and baryo-chemical potential at chemical freeze- 
out. We compare our results with the values obtained in other studies and also 
investigate more technical aspects such as a potential bias in the fits when fitting 
particle ratios or yields. Using parametrizations of the temperature and baryonic 
chemical potential as a function of energy, we compare the model calculations with 
data for a large variety of hadron yield ratios. We provide quantitative predictions 
for experiments at LHC energy, as well as for the low RHIC energy of 62.4 GeV. The 
relation of the determined parameters with the QCD phase boundary is discussed. 
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1 Introduction 



The success of the statistical (thermal) model [1] in describing the ratios of hadron yields 
produced in nucleus-nucleus collisions is remarkable. The thermal model was initially used 
for the AGS and SPS data [2] and was subsequently employed to describe data at SIS 
[3,4], SPS [5] and more recently at RHIC [6,7,8,9]. An analysis of the energy dependence 
of the thermal parameters extracted from fits of the experimental data, temperature (T) 
and baryo-chemical potential (/if,), established the "line of chemical freeze-out" [10]. These 
data were subsequently interpreted in terms of an universal condition for chemical freeze- 
out [11]. Remarkably, it appears that at the top SPS energy the (T,fib) values reach the 
phase boundary between the hadronic world and the quark-gluon plasma (QGP) [1,10], 
as calculated solving Quantum Chromo-Dynamics (QCD) on the lattice [12,13,14]. In this 
context, it has been argued that the QGP itself and the "deus ex machina" of phase space 
filling during hadronization are playing the crucial roles in achieving thermalization in 
high energy nucleus-nucleus collisions [15]. More recently it was demonstrated that pass- 
ing through the phase transition leads to multiparticle scattering of Goldstone bosons 



which drives even (multi)strange baryons rapidly into equihbration [16], providing a nat- 
ural explanation for the observation that the chemical freeze-out line reaches the phase 
boundary for small values of fib- The situation is less well understood for fit, >400 MeV 
and needs further investigation. 

Despite its simplicity, the thermal model allows to extract essential properties of the hot 
and dense fireball produced in high energy nucleus-nucleus collisions at a given stage 
in its evolution, namely when the inelastic collisions cease (chemical freeze-out). This 
makes the model a unique tool in the attempt to quantify from the experimental side 
the features of the phase diagram of hadronic matter [17,18,19]. Recent analyses over a 
broad energy range [20,21] have contributed in the on-going efforts to understand the 
chemical freeze-out criteria [1,10,11,22,23,24] in the phase diagram. It has been pointed 
out within a thermal model analysis that scanning the energy one encounters a transition 
from baryon- to meson-dominated freeze-out [25], with its associated fingerprint on the 
characteristics of hadron yields, also evidenced earlier [26]. 



Fig. 1. The energy de- 
pendence of experimen- 
tal hadron yields at 
mid-rapidity for various 
species produced in central 
nucleus-nucleus collisions. 
The energy regimes for vari- 
ous accelerators are marked. 
Note that, for SPS energies, 
there are two independent 
measurements available for 
the A hyperon yields. For 
more details see text (also in 
Section 3). 
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It is important to emphasize that experimental data on hadron yields are now available 
over a broad collision energy range and for a large sample of hadron species. A compi- 
lation of measurements of yields at mid-rapidity for the most abundant hadron species 
is shown in Fig. 1 for central nucleus-nucleus (Au-Au or Pb-Pb) collisions. As the cen- 
trality selection differs between various measurements, we have scaled the data at the 
same number of participating nucleons, Npart='i50. The contribution from feeding due to 
weak decays has been removed whenever the case (see Section 3). The main properties 
concerning the chemical composition of the fireball (at mid-rapidity) can be derived from 
the yields plotted in Fig. 1 without the need of any model. At lower energies {^/snn <5 
GeV), measured at Brookhaven's Alternating Gradient Synchrotron (AGS), the fireball 
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is dominated by the incoming nuclcons. while the yield of the produced pions, with a 
strong energy dependence, dominates at larger energies. The importance of isospin at 
lower energies is reflected in the different yields of 7r+ and n~. The decreasing yield of 
protons points to an increasing transparency of the incoming nuclei as a function of en- 
ergy. Beyond ^/snn —100 GeV, the newly produced protons become dominant. The yield 
of strange hadrons shows a sharp rise at AGS energies, with characteristic features for 
various species, determined by their quark content. The yields of and A (both with 
only the strange quark newly produced) are larger compared to K~ , which has two newly 
produced quarks. The remarkable similarity of the yields of K'^ and A, despite their 
large mass difference, is determined chieffy by their (anti)strange quark content, leading 
to their associated production, and less by the abundant presence of light quarks from 
stopped incoming nucleons. The yield of antiprotons and antihyperons (containing three 
newly produced quarks) is very similar and with a strong energy dependence (onset of 
production) at energies of CERN's Super Proton Synchrotron (SPS). At the Relativistic 
Heavy Ion CoUider (RHIC), due to a rather small net baryon content of the fireball, these 
differences almost disappear. The data are further discussed in Section 3 in comparison 
with thermal model calculations. We notice here that, in general, there is good consistency 
between overlapping data sets measured by different experiments. One exception are the 
hyperon yields at SPS energies, illustrated for A hyperons in Fig. 1. 

In view of these considerations it is important to assess with precision potential uncertain- 
ties in the extraction of thermal parameters with a consistent, 2nd generation, analysis 
of all available experimental data over a broad range of collision energies (y^J/viv=2. 7-200 
GeV) . This includes a comparison of the results from fits of hadron yield ratios and from 
absolute yields and a discussion on the infiuence of systematic errors on the fit parameters. 
The next section contains a brief description of the model. In Section 3 we perform fits 
of model calculations to the data. We discuss the energy dependence of the resulting T 
and Hb parameters, along with a critical assessment of various sources of uncertainties, in 
Section 4. Section 5 is devoted to the comparison of model and data in terms of excitation 
functions of various hadron yield ratios, which we extend up to the LHG energy. The 
relation of the extracted parameters with the predicted QGD phase boundary, calculated 
solving QGD on the lattice [12,13,14], is discussed in Section 6. The Appendix contains a 
discussion of more technical aspects in the fit procedures. 



2 Model description 

We restrict ourselves here to the basic features and essential results of the statistical model 
approach. A complete survey of the assumptions and results, as well as of the relevant 
references, is available in ref. [1]. 

The basic quantity required to compute the thermal composition of hadron yields mea- 
sured in heavy ion collisions is the partition function Z(T,V). In the grand canonical 
(GG) ensemble, the partition function for species i is {h — c — 1): 









from which the density is then calculated according to: 



AT TdlnZi gi 7 p^dp . . 



where gi = (2 Jj+1) is the spin degeneracy factor, T is the temperature and E,i = yp^ + mf 
is the total energy. The (+) sign is for fermions and (-) is for bosons. For hadron i of 
baryon number B^, third component of the isospin I^i^ strangeness S^, and charmness 
Cj, the chemical potential is /Xj = /ift-Bj + Hi^hi + IJ^sSi + i^cCi- The chemical potentials 
related to baryon number (/if,), isospin (///g), strangeness {^s) a-nd charm {yLc) ensure 
the conservation (on average) of the respective quantum numbers: i) baryon number: 
V^iTiiBi = Nb', ii) isospin: VJ2i''T'ihi = -^3°*; hi) strangeness: VJ2i'niiSi = 0; iv) charm: 
V J2i ^iCi = 0. The (net) baryon number Nb and the total isospin /g"* of the system 
are input values which need to be specified according to the colliding nuclei studied. 
The degree of stopping of the colliding nuclei, which is energy dependent and cannot be 
precisely determined experimentally, brings some uncertainty in the choice of Nb and Jg"*. 
In our case, as we study central collisions of heavy nuclei (Au or Pb), but focus on data at 
mid-rapidity, we have chosen Nb=200 and I^°^=-20. The sensitivity of the hadron ratios, 
which are the calculated "observables" to be compared with the experimental data, on 
Nb and /g"* is rather small. Taking into account the conservation laws i)-iv), T and the 
baryo-chemical potential Hb are the only parameters of the model, which will be obtained 
from fits to experimental data. When fitting hadron yields rather than ratios of yields, 
the fireball volume appears as additional parameter. We discuss its physical significance 
in Section 4 and in the Appendix. 

We want to emphasize that our model does not contain any strangeness suppression factor, 
75, as used e.g. in ref. [20,21]. Whenever 75 is used in thermal model calculations^, 
it is meant to account for non-equilibration in the strangeness sector. Departure from 
equilibrium is expected for elementary collisions or for peripheral nucleus-nucleus collisions 
(as well as for hght nuclei). The introduction of a 75 factor is an attempt to model this 
situation (see e.g. [27,20]). It is established that 75 ~ 1 for central collisions at RHIC 
[8,9]. At the lower end of the energy range, namely at SIS energies (up to 2 AGeV), there 
are indications [3] that 75 is not needed either. As we investigate central Au-Au or Pb-Pb 
collisions we prefer to stay within the thermodynamically well defined equilibrium model. 
Nevertheless we investigate, in the Appendix, the infiuence on our results as well as the 
statistical significance of the introduction of 75 into the fit procedure. This leads to the 
conclusion that our main results, i.e. the determination of T and /if, from fits of ratios of 
yields at mid-rapidity, are little changed by this additional fit parameter. 

The following hadrons are included in the calculations: i) mesons: non-strange (37), 
strange (28), charm (15), bottom (16); ii) baryons: non-strange (30), strange (33), charm 
(10); iii) "composites" (nuclei up to *^He and iir~-clusters [28], 18). The corresponding 
anti-particles are of course also included. Their characteristics, including a rather com- 
plete set of decay channels, (all strong decays and a suitable fraction of the weak decays 
contribution matching the experimental conditions) are implemented according to the 



We recall that the usage of 75 imphes that the thermal density of any given hadron carrying 
strangeness has a " suppression" factor 75 for every strange or antistrange quark. 



most recent PDG compilation [29]. Studies of QCD thermodynamics using lattice results 
have demonstrated that the partition function of a hadron and resonance gas model, which 
is a sum over all individual partition functions (Eq. 1), describes the lattice QCD data 
after suitably rescaling hadron masses to match lattice QCD [30]. 

We have used vacuum masses for all hadrons. The effect of hadron mass (and width) 
modification was analyzed in the context of the thermal model for the top SPS energy 
in ref. [31]. The conclusion is that at SPS a 10% reduction of masses leads to a slightly 
better thermal fit, with a corresponding 10% reduction of the resulting T and /Xb, A 
thermal analysis of the RHIC data using hadron masses derived from chiral models [32] 
found no preference for any scenario of in-medium masses compared to the non- interacting 
gas model. Note that at AGS (Si+Au) the nominal masses are preferred [1]. It was argued 
[33,34] that changes of masses in-medium in connection to chiral symmetry restoration 
can lead to larger values of the chemical potentials at lower energies. 

The finite widths of resonances are taken into account in the density calculation by an ad- 
ditional integration, over the particle mass, with a Breit-Wigner distribution as a weight: 

oo oo ry , 

^gi 3_ ( . f pMp 

NbwI i + rf/4 ■ exp[(£;r - ^Ji,)/T\ ±1' 



where rrii is the nominal mass and Fj is the width of particle i. The energy is then 
calculated for every value of m in the integration step, = \J-p^ + m^. Here, New is the 
normalization of the Breit-Wigner distribution and Mq is the threshold for the dominant 
decay channel. As this procedure significantly increases the computing time, for the fits it 
was used only for the AGS energies (see below), where in particular the widths of nucleon 
resonances arc expected to play a significant role [35] . At higher energies the width of the 
resonances plays little role in the determination of the thermal parameters. 

The interaction of hadrons and resonances is usually included by implementing a hard core 
repulsion of Van dcr Waals-type via an excluded volume correction. This is implemented 
in an iterative procedure according to ref. [36,37]: 

pexd.^j.^ /,) = p^'^-5-(T, /i); /i = // - KisenP'"^ (r, /x) (4) 



where V^igen is calculated for a radius of 0.3 fm, considered identical for all hadrons [5]. 
This correction infiuences the particle densities, but has httle (<5%) effect on ratios. The 

effect on absolute densities is very significant, however: at top SPS energies and above 
the excluded volume procedure leads to a reduction in densities by about a factor of 2, 
leading to a corrspondingly larger volume at chemical freeze-out. 

The grand canonical ensemble is the simplest realization of a statistical approach and 
is suited for systems with a large number of produced hadrons. However, for small sys- 
tems (or peripheral nucleus-nucleus collisions) and for low energies in case of strangeness 
production, a canonical ensemble (C) treatment is mandatory [38]. It leads to a phase 
space reduction for particle production (so-called "canonical suppression"). It has been 
shown [39,1] that the density of particle i with strangeness S calculated in the canonical 



approach, nf , is with a good approximation related to the grand canonical value, nf , 
as: rap = nf'^/Fg, with Fs = Io{x)/ Is{x). The argument of the Bessel function of order 
5" is the total yield of strange and antistrange hadrons. The suppression factor is shown 
in Fig. 2 for a canonical volume V(7=1000 fm^, which we employ in our calculations (see 
next section). In this case, the canonical suppression is negligible for all strange hadron 
species already for the highest AGS energy {^/snn —5 GeV) but is sizeable for the lower 
energies. Whenever the canonical suppression is needed, Vc, becomes another parameter 
of the model and needs to be evaluated carefully (see below). 
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Fig. 2. The energy dependence of the canonical suppression factor for strangeness, calculated 
for the canonical volume Vc'=1000 fm^. 

We mention that several similar codes are now available for thermal model calculations: 
SHARE [40], THERMUS [41], Therminator [42]. 



3 Fits to experimental data 

To extract the parameters of the model we perform fits of the experimental data with 
model calculations. In order to employ a minimal number of model parameters, namely 
T and Hb, the temperature and the baryonic chemical potential at chemical freeze-out 
(ceasing of inelastic collisions), a common practice is to fit hadron yield ratios [1,7,8,9]. 
It is important to recognize that one is faced with a choice when selecting — 1 sta- 
tistically independent ratios from the total number A^(A^ — l)/2 of ratios which can be 
constructed from A^ experimentally measured hadron yields. Although the number of in- 
dependent hadron yields is typically rather large (A^ >10) one may argue that fitting a 
particular choice of (statistically independent) particle ratios may introduce a bias for the 
extracted T and fib whenever a given experimental yield used in several ratios is subject 
to a systematic error . To minimize such a possible bias, we use ratios constructed from 



all measured hadron species avoiding, as much as possible, the repetition of any particular 
yield. Another possibility is to fit the measured hadron yields [20,21], which implies an 
extra parameter, the fireball volume. This method is free of the bias discussed above, but 
is subject of another bias, arising, for instance, whenever two measured yields (i.e. of a 
particle and its antiparticle) have a similar systematic error, which is cancelled in the 
ratio. Also, ratios typically have a smaller systematic error than absolute yields. Based 
on this, our main analysis focusses on ratios of hadron yields. As mentioned in Section 4 
and demonstrated in detail in the Appendix, the difference of fit parameters resulting 
from the two different fit procedures is always smaller than the uncertainty in the overall 
method, given the current status of measured yields and corresponding errors. 

We briefiy discuss below the fit procedures. The best fit is obtained by minimizing the 
distribution of x^. To get an estimate of the systematic error of the fits, we additionally 
consider the quadratic deviation, S^. The two quantities are defined as: 

( pexp _ r)therm\2 ( U^^P — ntherm\2 

2 _ \ " ) r2 _ \ ^ ) /r\ 



where R^^^ is the measured value of either the yield or the ratio of hadron yields ^ with 
its uncertainty o"j and i?*^^'"™ is the value from the model calculations. For the experi- 
mental errors we have quadratically added the statistical and the systematic errors of the 
measured yields. When the systematic errors of the measurements were not available we 
have assigned values of 10%, unless otherwise stated. The sum runs over the number of 
hadron ratios available experimentally. We focus on central collisions of Au or Pb nuclei. 
Ratios are calculated from the measured hadron yields. Whenever needed, we scale the 
measured yields with the number of participants [43] to account for different centrality 
classes. 

We focus on mid-rapidity data (dN/dy), for which the bulk of the published hadron 
yields is available, but we also analyze data integrated over the full solid angle whenever 
available. While at low energies the consideration of hadron yields integrated over An is 
optimal, this changes as, with increasing beam energy, the nuclei become transparent and, 
besides a central fireball, fragmentation regions develop. Once the stopping is not com- 
plete, necessarily the baryo-chemical potential depends on rapidity and a consideration of 
An data with one value of fij, would not be appropriate. In that case the statistical model 
should be compared to data over a central region near mid-rapidity. One may consider 
in separate analyses more forward or backward fireballs with different parameters. The 
possible presence of a finite net strangeness at mid-rapidity, as suggested by a recent 
transport model calculation [44], could pose a special difficulty; this is not accounted for 
in our model. However, the effect is expected to be small enough to not cause a systematic 
bias of the extracted thermal parameters. 

An important aspect in the comparison of calculations with measurements is the contri- 
bution of feed-down from weak decays, mainly in the yields of pions and ( ant i) protons. 
In our model, the fraction of those hadrons originating from weak decays is adjusted. In 
general, this contribution is subtracted in the experiments, either explicitly (via simu- 



2 We denote a yield ratio of two particles by the ratio of their respective symbols. 



lations) or implicitly (due to specific reconstruction methods). We have tried to model 
in detail the contributions from weak decays, following the information in the relevant 
experimental papers. In the Appendix we demonstrate the importance to determine as 
precisely as possible the size of this contribution. 



3.1 AGS 



At the AGS, collisions of Au nuclei at beam kinetic energies of 2 to 10.7 AGeV, corre- 
sponding to ys/v7v= 2. 70-4.85 GeV were studied. For the beam energies of 2, 4, 6, and 8 
AGeV, the yields of protons [45,46], pions [47,48] and kaons [48] are available (there are 
no K~ data at 2 AGeV). The A yields are only available integrated over Air [49] and for 
the energy of 6 AGeV the integrated S~ yield was also measured [50]. In these cases we 
have derived the mid-rapidity values assuming the same ratio between An and dN/dy as 
for protons [45]. 
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Fig. 3. Measured hadron yield ratios compared to the thermal model calculations for the AGS 
beam energies of 2, 4, 6, and 8 AGeV. The symbols are data, the lines are model calculations 
corresponding to the best fit (T, fib) values of (64,760), (78,670), (86,615), and (93,580) MeV, 
respectively. For the beam energy of 6 AGeV the ratios of 47r yields are plotted as diamonds. 

The comparison of the measured and calculated ratios for the best fit (x^ minimization) 
is presented in Fig. 3 for the four energies. Here and in the following plots we show both 



the K^/ti^ and the K~ /ti~ ratios, but in the fit either one or the other is included 
(as only one is statistically independent from the ratios tt~ /tt^ and K~/K~^), depending 
on the minimum value of x^- The resulting values of T and /ij, are between 64 and 93 
MeV and 760 and 580 MeV, respectively. The model describes the data very well; the 
values per number of degrees of freedom, x^/^df are 1.0/2, 0.43/3, 1.15/4, 1.14/3 for 2, 
4, 6, and 8 AGeV, respectively. Although this could be accidental, as the number of data 
points in the fits is very small, these low values could imply that the systematic errors 
of the measurements are actually overestimated. We mention that if we do not include 
in the calculations the finite widths of resonances the corresponding values of l^df are 
slightly worse, but the values of T and /i;, change only marginally. For both cases, the 
X^ distributions are narrow around the minimum and there is no correlation between T 
and as is observed for higher energies (see below). In some cases, close- by minima are 
identified and we have included this feature into the errors of the parameters. 
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Fig. 4. The contours (in steps of 1 from the minimum, marked by the dot) and yield ratios 
with best fit at mid-rapidity for the top AGS beam energy of 10.7 AGeV. The diamonds represent 
ratios of yields integrated over Att. The ratio K~ /'k~ was not included in the fit. Note the scaUng 
factor of 100 for the ratios p/p and A/A. 

At the top AGS beam kinetic energy of 10.7 AGeV (y^J/vw"=4.85 GeV) there is a large 
set of hadron yields experimentally available: p [51,52,46], d [52], vr [51,48], K [48,53], A 
[54], (j) [55], A [56]. The thermal fit to the data yields T = 124±3 MeV and /i^ = 537 ± 10 
MeV, corresponding to a minimum of x^/Ad/=6.9/7 (the errors correspond to 1 a). The 
model fits the data very well, as evident from the narrow x^ distribution and from the 
comparison of the ratios shown in Fig. 4. Not included in the fit is the ratio i^'~/7r~, which 
is not statistically independent from the ratios already included in the fit. The ratio A/A 
deviates from the model calculations substantially (A "anomaly" [56]). The experimental 
A/A is a factor of 2 higher than the model value and as such has a significant contribution 
to the x^ (if A/A is not included in the fit, the resulting values are: T = 123 ± 3 MeV, 
/if, = 538 ± 9 MeV, xVA^d/=4.0/6). 



The fit without the ratios c//p, p/p, A/A, and (f)/K^ gives T = 108 ± 9, /it = 555 ± 18 




Fig. 5. The canonical volume dependence of temperature extracted from thermal fits for the 
AGS beam kinetic energies of 2, 4, 6, 8, and 10.7 AGeV (in order of increasing temperature). 
Left panel is for minimizing x^, right panel is for (5^. The corresponding /Ndf and 5"^ /N^f are 
also shown. 



MeV, with x^/^d/=l-3/3. We note that, in this case, T and are not anymore correlated 
(as seen in Fig. 4), but rather loosely anticorrelated. We have performed this fit to test a 
possible bias in the extracted thermal parameters at the lower energies due to the smaller 
number of ratios experimentally available. We estimate in this way that the extracted 
temperature can be biased by up to 15%, while for /i^ the bias can be up to 3%. These 
systematic (upper limits) errors are used in the following. 

To investigate the sensitivity of the extracted {T,fih) values on the canonical volume Vc 
(which enters via the canonical suppression factor), we have performed the fits for all 
AGS energies for a range of Vc values. The results are presented in Fig. 5, where we show 
the dependence of T on the canonical volume parameter Vc for two cases: minimizing 
and 5^. Both cases result in very similar freeze-out temperatures. As expected, T slightly 
grows with decreasing Vc- The effect is gradually reduced towards higher energies and 
totally vanishes at the top AGS energy. Based on the dependence of /N^f and 5"^ /N^f 
on Vc (also shown in Fig. 5) one cannot constrain the canonical volume parameter. As a 
consequence, and since the canonical volume dependence is a rather small effect, we have 
chosen the value Vc'=1000 fm^. We mention that is unaffected by the canonical volume 
choice. 



3.2 SPS 



At the CERN SPS Pb-Pb collisions are measured in fixed target experiments, with beam 
momenta (f^energies) of 20 to 158 AGeV {^/snn =6. 27-17. 3 GeV). For the time being only 
the data collected at 40, 80 and 158 AGeV are final. The mid-rapidity yields of pions, 
kaons and (anti)protons and deuterons [57,58,59,60] are available for these energies. The 
yields of the complete set of strange (anti)hyperons [61,62,63,64,65] and of [66] are 
published for the top energy and at 40 AGeV. Since the deuteron yield is measured over a 
reduced range in momentum, we have decided to not include the deuterons in the fitting 
procedure for all the SPS energies, although the ratio d/p will be shown in the following 
plots. Also shown in the plots are the ratios S/vr" and fi/vr", which are not included 
in the fits because they are not independent from the ratios S/A and respectively, 
which are used because they result in a better fit. 
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Fig. 6. Distribution of and hadron yield ratios with best fit at the SPS beam energy of 40 
AGeV using the NA49 and NA57 mid-rapidity data (the ratios K~ /'k~ , d/p^ H/tt" and Jl/vr" 
are not included in the fits). The NA49 Att data are plotted for comparison (diamonds). For the 
ratios VL/tt^ and O/H, the Vt yield includes both Vt and Vt. 

At the beam energy of 40 AGeV (i/SiVAr=8.76 GeV), a combined fit to all data gives 
r=156l3 MeV, /ife=403li4 MeV, with xV^d/=15.4/13 (from the quadratic deviation 
analysis, r=166 MeV, /ib=434 MeV, 5^=0.58). The distribution and the comparison 
of the measured and calculated ratios for the best fit are shown in Fig. 6. A fit of NA49 
data alone leads to T=166±2 MeV, /i6=438±7 MeV, with x^/Ndf=&-2/7 (T=162 MeV, 
/i6=414 MeV, 52=0.10). 

The x^ distribution and the comparison of the measured and calculated ratios for the 
beam energy of 80 AGeV (y^s/viv=12.3 GeV) are shown in Fig. 7. The best fit is achieved 
for T=154 MeV, /ib=298 MeV, with 1 a errors of 6 and 20 MeV, respectively (T=152 MeV, 
Hb=271 MeV, with minimum 5^=0.23). The quality of the fit is not good, as seen from 
the value of x^/A'd/=23.4/5. We have noticed that the outcome of the fit is sensitive to 
the inclusion of different (statistically equivalent) hadron ratios. For instance, if the ratios 
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Fig. 7. Distribution of and hadron yield ratios with best fit at the beam energy of 80 AGeV 
using the NA49 data (the ratios K~ /ir~ and d/p were not included in the fit). The 47r yields 
(represented with diamonds) are plotted for comparison. 



p/tt~ and A/vr" are considered instead of p/vr" and A/vr" the resulting temperature from 
minimization is lower by about 10 MeV, while the temperature from 5^ minimization 
is unchanged. This shows the effect of the errors on the outcome of the minimization 
and indicates that the errors of the data are not entirely consistent. We have incorporated 
these results into the systematic errors of the thermal parameters. 

For the top SPS energy of 158 AGeV (y^s5v7v=17.3 GeV), there are mid-rapidity measure- 
ments available from the experiments NA49 [58,59,61,62,63], NA57 [64,65] (hyperons and 
Kg) and NA44 [57] (vr, K, p). Note that the NA44 pion yields are not explicitly corrected 
for the contribution from weak decays, but this is estimated to be small if one takes into 
account the reconstruction procedure in the NA44 spectrometer [67]. We have assumed 
an efficiency of reconstruction of 40% for pions from Kg (considered the average yield 
of K~^ and K~ measured by NA44) and of 20% from A (considered the average of the 
NA57 and NA49 yields) decays. In this way, 5.9 tt^ and 7.8 7i~ (dN/dy) are estimated to 
originate from weak decays. 

We have used in the fits the ratios S/A and fl/E. This choice leads to better values of 
compared to the case when using the ratios S/vr" and VL/tt~ instead. The ratio d/p was 
not included in the fits. The measured and calculated ratios for the best fit are shown in 
Fig. 8. A fit to all the combined NA49, NA57 and NA44 data gives T = 160 ± 5 MeV and 
/ift = 240 ± 18 MeV, with /Nfif='^Q/22 (from the least square deviation minimization, 
T=172 MeV, /i6=243 MeV, 52=0.86). For instance, if the ratios S/tt" and VL/tt- are 
considered instead of S/A and VL/'E, T = 155 ± 2 MeV and /if, = 224 ± 6 MeV, with 
xyNdf=Q3/22 (T=170 MeV, /ib=237 MeV, ^2=1.03). 

To investigate the contribution to these large values of x'^/^df, we have considered sepa- 
ratly two sets of hadron ratios: 
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Fig. 8. The ratios of hadron yields at mid-rapidity with best fits at the top SPS energy of 158 
AGeV. For the ratios of hyperons and Kg from NA57 relative to pions, we have used the tt~ yields 
measured by NA44. The 4tt ratios with NA49 data (diamonds) are plotted for comparison. The 
full lines are for the combined fit, the dashed ones for the NA44-I-NA57 data (T=180, /ifc=268 
MeV) and the dotted ones for the NA49 data (r=150, /ife=266 MeV). The ratios K'/tt", d/p, 
S/tt" and O/vr" are not included in the fit. 




(MeV) (MeV) (MeV) 



Fig. 9. The distributions at the top SPS energy of 158 AGeV for three cases (see text). 

i) the combined NA57 and NA44 data^ give T = 179 ± 7.5 MeV and /Xb = 267 ±26 MeV. 
With x'^/Ndf=8.7/10, the fit quality is good. 

^ In combining data from different experiments we neglect the possible difference in centralities 
at the same quoted fraction of the inelastic cross section which may arise from different centrality 
measures. 



ii) using the NA49 data wc extract T = 150 ± 4.5 MeV, = 226 ± 15 MeV, with 
/Nof— 27.8/10. The fit quahty is much worse compared to the NA57 and NA44 data 
set. 

The contours for the combined fit and for the two separate data sets are shown in 
Fig. 9. The outcome of the fits is summarized in Table 1. While fii, is compatible within 
the errors for all the three cases, the value of the temperature differs significantly, well 
beyond the errors, between the NA49 and the NA44+NA57 data sets. This discrepancy 
imposes large systematic errors for the extracted thermal parameters (see next section). 
Remarkably, the values extracted from the 5^ minimization are very similar for all three 
cases, T ~170 MeV, fi^ ^240 MeV. This indicates that the outcome of the combined fit is 
driven by the NA49 data and in particular by their small errors. Although unsatisfactory 
thermal fits may signal a possible change in physics in connection to the critical (end)point 
of a first order phase transition from hadrons to QGP [68] , the experimental situation just 
discussed allows unfortunately no conclusion on this issue. The apparent inconsistency of 
the results obtained using two independent data sets is not alleviated if one uses in the 
fit hadron yields instead of ratios. The results for this case (presented in Table 4 in the 
Appendix) show a similar pattern as in case of the fits using ratios. We note that poor 
quality of the fits at SPS energies characterizes other recent analyses within the thermal 
model [20,69]. 

Table 1 

Summary of the results of the fits using mid-rapidity data at top SPS energy (y^S7VJV=17.3 GeV). 
The first three and next three columns show the results of the (with their 1 a error) and 5^ 
minimization, respectively. 



data set 


T (MeV) 


/U6 (MeV) 




T (MeV) 


^ib (MeV) 


6^ 


NA44+NA57 


179±7.5 


267±26 


8.7/10 


174 


243 


0.15 


NA49 


150±4.5 


226±15 


27.8/10 


168 


240 


0.66 


combined 


160±5 


240±18 


56/22 


172 


243 


0.86 



3.3 RHIC 



A complete set of final data for Au-Au collisions at RHIC is available at y^J/vw"=130 GeV. 
We use for our analysis the yields of pions, kaons, (anti)protons [70] and A, A [71] measured 
by PHENIX, while for the other hadron species, S, Q [72,73], [74], K* [75], the data are 
from STAR. Most of the available data are not corrected for feed-down from weak decays, 
but these contributions are properly taken into account in the model calculations. One 
exception are the pions, for which the weak decay feed-down contribution is not clearly 
specified by the experiment (PHENIX) .We have considered three cases: 

i) Assuming that all the pions from feed-down are included in the data, the values of the 
thermal parameters are T — 170.5 ± 6 MeV and = 38 ± 15 MeV, with a minimum 
X'^/Ndf value of 3.6/11. 

ii) At the other extreme, if one considers that the pion data do not contain any contribution 
from weak decays, the resulting values are T — 162.5 ± 5.5 MeV, /Xf, = 35 ± 11 MeV, 
xV^ci/=4.6/ll. 

ni) Assuming that the data contain 30% of the pions from weak decays, T — 165.5 ± 5.5 
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Fig. 10. The distribution and hadron yield ratios with best fit at y/sNN=130 GeV. Here, the 
Q yield includes both $7 and 

MeV, /ife = 38 ± 11 MeV, /Ndf=4:.l/ll. 

Although the fit is good for all cases (and even if the minimal is obtained for case i), we 
consider the intermediate case iii) as the most likely situation, as it could be an implicit 
result of the reconstruction in the experiment. The outcome of the fit is shown in Fig. 10. 
The very good fit is also apparent in the comparison of the hadron ratios: essentially all 
the experimental ratios are well reproduced by the model, including those involving and 
K* resonances. 

At y^Sjvjv=200 GeV only a limited set of yields of identified hadrons is available to date. 
The yields of tt^, K^, p, and p are published by PHENIX [76], STAR [77], and BRAHMS 
[78] and the values agree within the quoted systematic errors. In addition, available are 
the yields of [79,80], ir(892)* [81], d and d [82]. We use the hadron ratios with the 
corresponding errors whenever provided by the experiments [76,79,80,81] or otherwise 
calculate the ratios using the published yields quoted above. Further ratios are available 
as preliminary data on strange hyperon ratios [83], A++/p [84], p/vr~, A/tt", S/tt^, Q/tt~ , 
and A* /A [85]. Unless specified, our fits do not include the ratios involving resonances 
(A++/P, K{892)*, and A*). 

We consider the following cases for the fits: 

i) a combined fit of all available data, with the exception of strongly decaying resonances: 
T = 155 ± 2 MeV, /Xb = 26 ± 5 MeV, xV^d/ = 34.1/23 (with 5^ minimization: T=164 
MeV, /i;,=24 MeV, 5^=0.40). If we include in the fit the three ratios involving resonances 
{K*/K~, A*/A, and A++/p), the results are the same within the errors, T = 155 ± 2 
MeV, /ib = 25 ± 5 MeV, but with a worse xV^df = 41.8/26 (T=162 MeV, /ife=22 MeV, 
5^=0.82). 

ii) as i), but excluding from the fit the ratios p/tt~ and (p/K' from PHENIX. The result- 
ing parameters are: T = 160.5 ± 2 MeV, /i^ = 20 ±4 MeV, with xV^df =16.0/21 (T=166 
MeV, /i6=26 MeV, 5^=0.19). Including the resonances the outcome of the fit is identical, 
T = 160 ± 2 MeV, /Xb = 20 ± 4 MeV, with a reasonable xV^df = 25.2/24 (T=164 MeV, 
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Fig. 11. Hadron yield ratios with best fits at y^sjv7v=200 GeV (see text). The last three ratios, 
involving resonances, were not included in the fits. 

fib=20 MeV, 5^=0.61). 
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Fig. 12. Distributions of at y/sNN=200 GeV for the two cases in Fig. 11: a combined fit of 
all available data (left panel) and excluding the ratios p/iT~ and (jj/K^ from PHENIX (right 
panel). 



Obviously, the outcome of the fit (lower T) for case i) is determined by two contributions: 



a) the 6% systematic error assigned by PHENIX to the ratio p/i^' ^ which is much smaller 
than the values 14% and 18% which we calculate from the systematic errors of the absolute 
yields in case of BRAHMS and STAR data, respectively; b) the ratio (t)/K~ from PHENIX 
[80] , which is much lower compared to the value from STAR [79] . Using all data from STAR 
(without resonances) and the d/p and d/p ratios from PHENIX gives T = 160.5 ±2 MeV, 
/ift = 15 ± 7 MeV, with xV^df = 11-6/13 (from 5^ minimization: T=166 MeV, /ib=22 
MeV, 5^=0.17). Taking into account that, for all cases, the 5^ minimization results lead to 
T >160 MeV, and despite the fact that the x^/^df is reasonable even for the global fit, we 
consider the fit without the ratios p/vr" and 4>/K~ from PHENIX (case ii) as the nominal 
case for ys/vjv=200 GeV. The systematic errors are then estimated as the differences to 
the case when all measured hadrons (but no resonances) are included in the fit and to the 
values from the 6"^ minimization. 

The experimental ratios and the best fit calculations for the two cases, i) and ii) are shown 
in Fig. 11. Their respective distributions are shown in Fig. 12. We note that, comparing 
the resonances, only the A* ratio is clearly deviating from the model fits in our case, at 
variance with other claims in this respect [86] . Despite this deviation, which needs further 
support from the experimental side, our results on resonances do not point to a sizeable 
hadronic rescattering effect after chemical freeze-out [86]. 
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Fig. 13. The dependence of various hadron ratios on T and fib for the range of parameters 
relevant for the RHIC energies. Note that the ratios H^/vr^ and /ir^ are scaled up by a 
factor 10. 



In Fig. 13 we show the sensitivity of various hadron ratios on T and fib for the range 



of values relevant for the RHIC energies [6]. The importance of particular ratios in con- 
straining T and/or yUf, is evident. The antihyperon to hyperon ratios as well as the ratios 
of heavy hadrons to pions are of special value in this respect. 

We note that thermal models have also been used to describe hadron production in e+e" 
and hadron-hadron collisions [38,27], leading to temperature parameters in the range 160- 
170 MeV. We mention that our model does not fit well the e+e^ data. Our fits to the 
data at y^SiV7v=29 GeV compiled in ref. [27] were performed constraining Hb=0 and with 
three free parameters: the temperature, T, the strangeness suppression factor, 75, and the 
volume, V (used to calculate the canonical suppression factor). The fits yield T — 154 ±3 
MeV, 75 = 0.74 ± 0.06 and V = A2 ± U fm^, with a xViVd/=56/ll. Obviously, the fit 
quality is very poor. For comparison, in the analysis of ref. [27] T = 159 ± 2.6 MeV, 
75 = 0.71 ± 0.05, x^/A'd/=29.3/12. Both T and /ih are compatible within the errors, but 
our fit shows a much worse quality. This may be due to the approximations concerning 
canonical strangeness suppression which are employed in the present approach. However, 
even with the appropriate canonical treatment as done in ref. [27], the quality of the fits 
(both in e"'"e~ and in hadron-hadron collisions) is not good, raising questions about the 
applicability of the statistical model for elementary collisions. 



4 Energy dependence of the thermal parameters 

Our results for the best fits are summarized in Table 2. The upper part in Table 2 shows 
the results of the fits using mid-rapidity data, the lower part is for the data integrated 
over 4:71. One notices that in general x^/^df is close to (and often below) unity, being 
significantly above unity only for the SPS energies. The fits 47r-integrated data result in 
comparable temperature values as for the mid-rapidity data in case of the AGS energies 
and sizably lower values for the SPS energies. For the higher SPS energies, the fits of Att 
data are characterized by even poorer values. We would like to stress that, although 
we have performed the fits of An data up to the top SPS energy (in order to compare to 
similar analyses in other works) , we consider the extracted thermal parameters in this case 
adequate only up to the energy of 40 AGeV. Also included in Table 2 are the results of the 
5^ minimization, which differ sometimes from those obtained from the minimization. 
The central values used in the following are those from the minimization. 

The " systematic" errors AT*^** and A^l^^^ shown in Table 2 are obtained by quadratically 
adding the 1 a error from the minimization with the difference between the values 
obtained by x"^ and 5^ minimization . Whenever redundant sets of data are available from 
different experiments, which is the case at top SPS and top RHIC energies, our "nominal" 
values are those resulting from the combined fit of all data. We have assigned as systematic 
error the quadratic sum of 1 cr fit errors and the deviation between the (7",//^) values of 
this case and those resulting from the fit of the data which results in the smallest value of 
(see previous section and Table 2). For the SPS energy of 80 AGeV {^/sJ^=12A GeV) 
the systematic errors are evaluated based on the sensitivity of the fit on the inclusion of 
different hadron ratios. As discussed in the previous section, for lower AGS beam energies 
of 2-8 AGeV, we have estimated 14% and 3% systematic errors for T and Hb, respectively, 
based on the sensitivity of the fit using ratios at top AGS energy. For the AGS energy 



Table 2 

Summary of the results of the thermal fits using ratios. The upper part shows the results obtained 
using mid-rapidity data, the lower part is for data integrated over 4tt. The first three and next 
three columns show the results of the (with their 1 a error) and minimization, respectively. 
The quantities AT*^** and A/x^^** (in MeV) are the assigned systematic errors (sec text). 
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of 6 AGeV, we have used the available An yields of pions [47], protons [45], A and 
[49]. We note that, at this energy, there are no K'^ and K~ yields measured in full phase 
space. Based on a comparison of the ratio K~ / at mid-rapidity and in 47r at top AGS 

energy [53], we have assumed that at 6 AGcV the ratio K~ / in Atx is 10% smaller 
than at mid-rapidity. At the top AGS energy, some of the 47r yields do not correspond to 
complete coverage [53]; also, the ratio p/p is obtained with a rather large extrapolation. 

In Fig. 14 we show the energy dependence of T and fi^ extracted from our thermal fits. The 
temperature T exhibits a sharp rise up to a/J/vw —7-8 GcV, while ^Uf, sharply decreases 
all the way up to RHIG energies. Our results are compared with results from other stud- 
ies [1,3,8,20,21,69], which exhibit similar trends. One notable exception are the results 
of Letessier and Rafelski [21]; in this case, the observed difference may arise from the 
usage in their work of seven free parameters, including, besides strangeness fugacity and 
suppression factor (75), a light quark occupancy factor (7^) and an isospin fugacity. The 
non-monotonic change in 7^ and 75 as a function of energy determines the temperature 
extracted in ref. [21]. To alleviate the poor quality of the fits at SPS energies, the model of 
Dumitru et al. [69], introducing an inhomogeneous freeze-out scenario, goes beyond other 
thermal models. 
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Fig. 14. The energy dependence of temperature and baryon chemical potential. Left panel: the 
results of the present work are compared to the values obtained in earlier studies (see text); 
right panel: the results of our fits using dN/dy data, both with ratios and yields. The lines are 
parametrizations for T and fib (see text). 

Not surprisingly, our results support those obtained earlier within the same model imple- 
mentation [1] . In detail, there are differences between our results and the rest of the other 
results. Our values for the analysis of the Att data are significantly different from those 
of Becattini et al. [20], obtained within a model employing the strangeness suppression 
factor. A remarkable agreement between our results and the analyses in ref. [6,8] is seen 
in case of the RHIC data, in particular at ^/sNN=130 GeV. The results at yj/v7v=200 
GeV [8,9] are also in agreement with our values. The higher T value at ^/S]\^N=200 GeV 
in ref. [6] is due to preliminary data. 

We have parametrized our results from the fits of mid-rapidity data (left panel in Fig. 14 
and Table 2) as a function of ^/snn (in GeV) with the following expressions: 



T[MeV] = Tu^ 1 



0.7 + (exp(y5^(GeV)) - 2.9)/1.5^ 



l + &v^(GeV)' 



(6) 
(7) 



where the parameters a = 1303±120 MeV and b = 0.286±0.049 GeV ^ are the results of 



a fit (x^/Aid/=0.48/8). Our fit parametrization is the one proposed in ref. [26], but with 
different parameters to better fit the newly obtained Hb values of the present analysis. In 
Eq. 6 we assume a "limiting" temperature Tum, which was obtained by fitting the five 
points for the highest energies (SPS and RHIC). The result of the fit is = 161 ± 4 
MeV, with /Ndf=0.3/3. 
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Fig. 15. Energy dependence of the volume for central collisions (A'part=350). The chemical 
freeze-out volume for one unit of rapidity, dV/dy, is compared to the kinetic freeze-out volume 
from HBT measurements, Vhbt [87]. Note that the scales are different for the two observables. 

We now briefly turn to another interesting parameter which is (either implicitly or ex- 
plicitly) determined in the course of thermal model analyses. The volume at chemical 
freeze-out (corresponding to a slice of one unit of rapidity, dV^/dy) is shown in Fig. 15 as 
a function of energy. The values extracted directly from the fits of particle yields (see Ap- 
pendix) are compared to the values obtained by dividing measured charged particle yields 
with calculated densities (based on the above parametrization of T and /i^; note that for 
the AGS energies of 2-8 AGeV the values of T corresponding to the upper limit of the 
systematic error were used instead). As expected, the two methods give identical results, 
with the exception of a small discrepancy for the lowest energies. The chemical freeze-out 
volume is compared to the kinetic freeze-out volume extracted from Hanbury Brown and 
Twiss (HBT) measurements, Vhbt [87]. We note here (see also the Appendix) that to 
determine a volume from thermal model analyses one needs to know absolute densities. 
This implies an explicit dependence of the volume on whether or not excluded volumes 
are implemented in the calculations, and, if included, on details of the implementation. 
We follow here the procedure developed in [5]. 

While the bias towards unphysically large volumes seen at the energies of 2-8 AGeV (see 



Appendix) is not completely understood, it is clear that the volume at kinetic freeze- 
out should not be exceeded. Note that, to relate quantitatively the magnitude of the 
rapidity density of the chemical freeze-out volume dV/ dy to the volume at kinetic freeze- 
out determined from HBT measurements, one needs to map rapidity onto space (see the 
discussion in [87]). Here we are mainly interested in the energy dependence of the two 
volumes. It appears that the volume at chemical freeze-out does exhibit a similar non- 
monotonic behavior as the volume at kinetic freeze-out, a remarkable result considering 
that the latter is determined by a completely different procedure. A minimum is observed 
for the top AGS energy, followed by a logarithmic increase as a function of y^iivw- 



5 The energy dependence of hadron ratios 

Using the parametrizations of the chemical freeze-out parameters derived from the fits 
of experimental data at mid-rapidity, as shown in Eq. 6 and 7, we calculate the energy 
dependence of various hadron ratios. No contribution from weak decay feed-down is in- 
cluded in the model calculations as the published data are mostly corrected in this respect. 
Whenever not the case, the measurements are scaled appropriately (see below). The fi- 
nite widths of resonances are taken into account in the model. This contributes about 
10% of pions at low energies and 4% at higher energies. Within the smoothing hypothesis 
implied by the parametrizations of Eq. 6 and 7, the model has interpolative and extrap- 
olative predictive power. In this sense, we can provide quantitative predictions for LHC 
energy, as well as for the low RHIC energy of 62.4 GeV, for which the experimental data 
are becoming available [88] . This also applies to the low energy SPS data [89] . 

In Fig. 16 we show the energy dependence of antiparticle to particle ratios. Calculations 
employing the global T(^S7VAr) and /Xb(y^iivjv) from Eq. 6 and 7 reproduce very well 
all these experimental ratios, including the preliminary values for ^/sNN=62A GeV [88]. 
The 7r~/7r"'", K~ /K~^, and p/p ratios illustrate the decreasing importance of isospin, the 
increasing importance of strangeness production and the dramatic evolution as a function 
of energy of the balance between the incoming and newly produced protons. The ratio 
p/p has the steepest dependence on energy. Prom the lowest measured value of 3-10~^ at 
AGS, it reaches a value of about 0.7 at RHIC and is predicted to reach 0.95 at LHC. The 
dominance of the initial nucleons at low energies changes to the newly created hadrons 
at RHIC and beyond. At LHC, it is expected that the fireball will consist practically 
completely of created hadrons (//(, ~1 MeV). 

The energy dependence of antihyperon/hyperon ratios (right panel in Fig. 16) follows a 
mass hierarchy: the saturation value of 1 is achieved the earlier the more massive the 
hyperon species. This is a fingerprint of the preferred abundance of hyperons over anti- 
hyperons containing 2 (A) or 1 (5) light valence quarks, which may be remnants from 
the incoming nucleons, while in case of Q, all valence quarks are newly produced. ^ Re- 
markably, as already seen in Fig. 8, there is a very good agreement between the NA49 
[61,62,63] and the NA57 [64] data. The ratio A/A from the NA45 (CERES) experiment 

^ The small structure of the excitation function at ^/sNN=10-20 GeV, most prominently seen 
in Fig. 16 for the ratio Cl/n, but also for S/S, is an artifact of our parametrization of fif,. 
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Fig. 16. The energy dependence of antiparticle/particle ratios. 

[90] is as well in agreement. The generally poor quality of thermal fits at SPS energies is 
not caused by ratios involving (multi) strange baryons. 

The comparison of the measured and calculated excitation functions for hadron abun- 
dances with respect to pions is shown in Fig. 17. For a consistent comparison to the other 
measurements and to the calculations, the measured pion yield at i/J/vw"=130 GeV [70] 
was scaled assuming that 30% of the weak decays contribution was contained in the data 
(see above). The ratios shown in Fig. 17 have a very different energy dependence com- 
pared to those in Fig. 16, reflecting the evolution of the fireball composition at freeze-out 
as a function of energy. The steep decrease of the p/vr^ ratio directly reflects the decrease 
as a function of energy of stopping of the incoming protons, implying a decrease of Hb- 
The increase of pion production also plays a role in this ratio. Beyond -s/s^n —100 GeV, 
the flattening is a consequence of the dominance of newly created baryons. The steep 
variation of the K~^/tt^ and K~/7t~ ratios at the lowest energies reflects the threshold 
for strangeness production, determined in the model by the steep increase of the temper- 
ature. The canonical suppression plays an important role as well. While the ratio /n^ 
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Fig. 17. The energy dependence of hadron yields relative to pions. In the right panel, note the 
scaling factors of 10 and 100 for the H^/vr^ and ri^/vr^ ratios, respectively. 

shows a monotonic increase with energy, followed by a saturation, essentially determined 
by the temperature (as both particles are newly created), the ratio /t:^ , which is much 
discussed as a signature for the onset of QGP [91,92,93,94], shows a characteristic broad 
peak around ^Jsnn —8 GeV. As the contains a u valence quark, which may come 
from the initial nucleons, its yield is the convolution of two competing contributions as 
a function of energy: i) the decreasing net light quark content and ii) the increasing pro- 
duction of quark-antiquark pairs. The peak in the /tt^ ratio occurs naturally in the 
thermal model [26,25], but is broader and has to be seen in the context of other strange 
hadron yields [25] (right panel of Fig 17). It appears that, at the SPS energies of 80 and 
158 AGeV, the yields of kaons and protons relative to pions are systematically below the 
thermal model predictions. We note that recent transport model calculations [93] grossly 
deviate from the measured /tt^ ratio, while earlier results reproduced the data well 
[95]. We emphasize the good agreement between the two data sets (NA49 [58] and NA44 
[57]) available at the top SPS energy. 



The energy dependence of the relative hyperon yields, A/vr", "E' /t[~ and VL~ /ti~ ^ shown in 
the right panel of Fig. 17, reveals the presence of characteristic peak structures, already 
noted in [26]. Their strength and location follows a mass hierarchy, recently discussed 
by Cleymans et al. [25]. The peaks are less pronounced and located at larger energies 
{y/SNN — 5, 10 and 20 GeV) the more massive the hyperon species. This results from 
an interplay between the baryochemical potential (presence of the light quarks from the 
initial nucleons) and temperature. The agreement between the model and the data is good 
at AGS and RHIC energies, while at SPS the situation is more complex. The two available 
sets of data, NA49 [61,62,63] and NA57 [64], are in a clear disagreement: generally, the 
NA49 data are below the NA57 data. As a consequence of our adopted fitting procedure 
(see above) the model prediction is between the two data sets. We notice that the discrep- 
ancy of the data is manifest in the absolute yields: for instance, dN/dy of A at top SPS 
(5% central collisions) for NA49 and NA57 are 11.0±1.6 and 18.5±2.2, respectively. The 
disagreement is well beyond the errors (statistical and systematic, added in quadrature) 
quoted by both experiments. 




Fig. 18. The energy dependence of d/p and (f)/K ratios. 

In Fig. 18 we show the excitation function of the ratios d/p and (f)/K~. The ratio d/p 
shows a monotonic decrease as a function of energy, well explained by the model. The 
ratio (()/ K~ shows a peak in the calculations, centered in this case at ^/snn —8 GeV. The 
model reproduces well the AGS measurement [55], while at SPS there is a disagreement. 
At RHIC the experimental situation is unclear; the model is in agreement with the STAR 
data. 

Given the accuracy of the description of hadron ratios presented in the previous section, 
it is not surprising that, based on a parametrized energy dependence of T and /ib, the 
model does reproduce the experimental values well over a broad energy range. However, 
systematic discrepancies between the model calculations and data are evident for the SPS 
data, in particular concerning the ratios involving strange to non-strange hadrons. 
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Fig. 19. Energy dependence of strangeness to entropy ratio (see text). 

A more global way to represent the ratio of strange to non-strange hadrons is the strangeness 
(ct) to entropy {S) ratio. Its excitation function is presented in Fig. 19. We adopt an 
experiment-oriented construction of the two quantities, which we consistently employ for 
the model as well, calculated from the yields at mid-rapidity as: 

a = 2 X + + 1.54 X (A + A), 5 = 1.5 x (7r+ + tt") + 2 x p (8) 

The strangeness has in principle to be complemented with the yields of 0, S and Q (and 
H and Q), but, since the measurements for these yields are scarce (as seen in Fig. 17 and 
18), we have chosen to leave them out for the strangeness count. The factor 2 multiplying 
the kaon yields takes into account K*^, while the factor 1.54 for A hyperons accounts for 
the contribution of and was deduced from the model calculations. ^ The factor 1.5 for 
the pion yields accounts for the 7r° yield, while in case of p yields the factor 2 is used to 
account for the produced protons. As expected from the individual particle ratios studied 
above, the ratio strangeness/entropy is well reproduced by the model, with the exception 
of the data at SPS, where the NA49 data exhibit a sharper peak than predicted by the 
model. This feature is not supported by the NA57 data, which are in agreement with the 
model. 

The results discussed above allow a rather confident prediction for hadron ratios in central 

^ The yield of T,^, which decays with 100% branching ratio into A7, is always included in the 
A yield. 



Table 3 

Thermal model prediction for hadron ratios at LHC (T=161 MeV, //(,=0.84 MeV). No contri- 
bution from weak decays is included. 
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Pb-Pb collisions at LHC energy [96], presented in Table 3. The antiparticle/particle ratios 
are all very close to unity, with the exception of p/p, reflecting the small, but nonzero, iij, 
value obtained from Eq. 7. Within the extrapolation scenario for the thermal parameters, 
the model predicts that the yields relative to pions are very similar to the values measured 
at the RHIC energy. 



6 The phase diagram of hadronic matter 

The values of T and obtained from our fits of the experimental data are shown in a 
phase diagram of hadronic and quark-gluon matter in Fig. 20. Full points are from fits 
of ratios of yields at midrapidity, open points are from fits of 47r ratios. An important 
observation about the phase diagram is that, for the 40 AGeV SPS energy and above, 
the thermal parameters agree with the phase boundary calculations from lattice QCD 
(LQCD) [12,13,14], implying that hadron yields are frozen at the phase boundary. The 
LQCD calculation [12] shown in Fig. 20 is for two light quarks (m, d) with realistic (close 
to physical) masses and a heavy strange quark. The critical temperature at ;Li5=0 from 
LQCD calculations is about 165 MeV [18] (and refs. therein), with a scale uncertainty of 
the order of 10 MeV [18,14] and with comparable systematic errors [18]. 

Also included in Fig. 20 are calculations of freeze-out curves for a hadron gas at constant 
energy density (e=500 McV/fm'^) and at constant total baryon density (n5=0.12 fm~^) 
[22]. The LQCD phase boundary calculated in [12] and shown in Fig. 20 apparently does 
not follow the expectation [13] (see also [97]) of a constant energy density dependence. 
However, depending on the method, the errors for the critical line from LQCD calculations 
can be large for finite [13]. 

The freeze-out points which are departing from the LQCD phase boundary are approx- 
imately described by the curve of a hadron gas at constant baryon density. An earlier 
proposed freeze-out criterion corresponds to an average energy per average number of 
hadrons of approximately 1 GeV [11]. It has also been argued recently that the freeze- 
out points can be described by a constant entropy density divided by [23]. This is 
a measure of the degrees of freedom; in a hadron gas with reahstic masses the effective 
number of degrees of freedom depends explicitly on temperature and is not constant. As 
Cleymans et al. [25] have recently pointed out, for realistic excluded volume corrections. 
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Fig. 20. The phase diagram 
of hadronic and quark-gluon 
matter in the T-^b plane. The 
experimental values for the 
chemical freeze-out are shown 
together with results of lattice 
QCD calculations [12]. The 
predicted [12] critical point is 
marked by the open trian- 
gle. Also included are calcula- 
tions of freeze-out curves for a 
hadron gas at constant energy 
density (e=500 MeV/fm^) and 
at constant total baryon den- 
sity (nf,=0.12 fm-3). The full 
triangle indicates the location 
of ground state nuclear matter 
(atomic nuclei). 
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this criterion grossly deviates from the freeze-out points. A comparison of various freeze- 
out criteria was recently done in ref. [24], which showed that all are identical except for 
large and small Hb values. All are smooth curves and consequently not consistent with the 
new results presented here, which exhibit a rather steep trend at intermediate fib values. 
However, the errors need to be improved before one can confidently rule out any (smooth) 
universal freeze-out criterion. An exciting possibility is that the rather abrupt turn-over 
in the freeze-out points near /ib=400 MeV is caused by the approach to the QCD phase 
boundary. 

The underlying assumption of the thermal model used to extract the (T,/ib) values is equi- 
librium at chemical freeze-out. A natural question then is how the equilibrium is achieved? 
The answer obviously cannot come from within the framework of the thermal model. Con- 
siderations about collisional rates and timescales of the hadronic fireball expansion imply 
that at SPS and RHIC the equilibrium cannot be established in the hadronic medium 
below the critical temperature Tc [15,16]. In a recent paper [16] many body collisions near 
Tc were investigated as a possible mechanism for the equilibration within the hadronic 
stage. There it is argued that because of the rapid density change near a phase transition 
such multi-particle collisions provide a natural explanation for the observation of chemical 
equilibration at RHIC energies and lead to T = to within an accuracy of a few MeV. 
While this argument is expected to be valid also for SPS energies, the situation at lower 
energies needs further consideration. 



The critical temperature determined from RHIC and SPS data assuming T ^ Tc coincides 
well with lattice estimates [12,14,18]. These arguments buttress the conclusion that the 



hadrochemical frcczc-out parameters probe experimentally in a unique manner the critical 
line of the QCD phase transition between hadrons and QGP [98]. Despite the rather large 
systematic errors of the extracted temperature for the SPS energies, our results imply that 
the phase boundary is reached for beam energies around 40 AGeV. The experimental data 
at 20 AGeV [89] will constrain further this energy range. 

Based on the LQCD results of Fodor and Katz [12] shown in Fig. 20, the experimental 
freeze-out points at SPS are located in the vicinity of the critical point. It was pointed out 
recently [19] that the existence of a critical point for Hb <500 MeV requires a fine tuning 
of the (light) quarks masses within 5%. However, it is important to recognize that serious 
open problems of LQCD [19] need to be solved before one could address quantitatively such 
a delicate possibility. Nevertheless, it is interesting to speculate whether the deviations 
from the thermal model (including rather poor-quality fits) which we have encountered for 
the SPS energies are a hint for the critical point. It is expected that, in the (broad) vicinity 
of the critical (end)point the thermal model would not work [68]. Thermal fits including 
fluctuations have been already performed for the top SPS energy [69]. Unfortunately, the 
present experimental situation, namely the level of disagreement between data, does not 
allow any firm conclusion on the interesting issue of the critical point. 



7 Summciry 

We have analyzed the experimental hadron yields ratios over a broad energy range 
(\/siV7v=2. 7-200 GeV) in comparison with thermal model calculations. The fits of the 
experimental data with the model provide the temperature and baryo-chemical poten- 
tial, at hadrochemical freeze-out. The quality of the fits is very good for most cases, 
providing support for the validity of the approach. From our analysis we have estab- 
lished parametrizations of T and iij, as a function of energy. This allows to assess and 
understand trends visible in the experimental data. In particular, we have discussed for 
the yields of strange and multi-strange hadron yields relative to pions the resulting non- 
monotonic energy dependence, which describes the main features observed in the data, 
but not the details. Furthermore, we have provided quantitative predictions for hadron 
ratios at y^SiV7v=62.4 GeV and at LHC energy, where first data are expected soon. The 
analysis presented here has also led to an updated version of the chemical freeze-out curve. 
The results buttress previous observations that, from -y/sNN ~ 10 GeV on, the chemical 
freeze-out points coincide with the phase boundary between hadrons and the QGP as 
predicted by solving QCD on the lattice. We have provided chemical freeze-out points, 
for the first time, at the 4 low AGS energies. In this energy range the freeze-out curve 
exhibits a peculiar and intriguing structure which needs further explanation. 



Appendix 

In this Appendix we discuss technical aspects related to the thermal model analysis. 
First, we concentrate on the technical aspects of fitting yields rather than ratios, and the 
demonstration that both methods produce equivalent results. Secondly, we briefly address 



the issue of the role of 75 in the thermal model fits. Finally, we comment on the issue of 
weak feeding. The absolute yields (and consequently the values for the volume) quoted 
here are for A'part=350. 



A. Analysis of particle yields 



The fit of mid-rapidity yields for the top ASG energy is shown in Fig. 21. The resulting 
parameters are: T=132 MeV, /ib=546 MeV, V=900 fml With xV^d/=10.3/7, the fit is 
good. The temperature is somewhat higher than the temperature obtained when fitting 
ratios. Excluding from the fit the yields of p, A, 0, and d one obtains: T=110 MeV, /i;,=550 
MeV, 1/=2620 fm^, xV^d/=l-2/3. The resulting change in T is 20%, larger than the 
difference derived from the fit of ratios (15%). This lower temperature is "compensated" by 
an unphysically large volume, leading to a good fit in this case too. The main contribution 
to this effect arises from the absence of p and A (when fitting without p and A: T=114 
MeV, /ife=551 MeV, 1^=2100 fm^, xViVd/=2.8/5). 
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Fig. 21. Hadron yields at mid-rapidity and the best fit for the top AGS energy. 

Our fits of (mid-rapidity) yields at lower AGS energies result in identical values for T 
and Hb as when using ratios. The extracted volume is much larger than for the top AGS 
energy. This originates in the bias of the temperature due to the limited set of particle 
yields available, as discussed above and as also observed when fitting ratios. Based on 
the above results at top AGS energy, we conclude that the extracted volume at the lower 
energies is overestimated by a factor of 2.9. 

The comparison of mid-rapidity yields and the best fit for the top SPS energy is shown 
in Fig. 22. The temperature is slightly higher than in case of ratios, but this is within 1-a 
errors. The discrepancy between different data sets noted in case of ratios is evident also 



Table 4 

Summary of the results of the fits using mid-rapidity yields at top SPS energy. 
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Fig. 22. Mid-rapidity yields and best fit (combined data) at the SPS beam energy of 158 AGeV 
[d was not included in the fit). 

for the yields. As for the ratios, this leads to contradictory results when fitting the data 
separately, see Table 4. The quality of the fits is for all SPS energies significantly poorer 
than in case of ratios. 

At the RHIC energy of i/J/vw"=130 GeV, Fig. 23, the resulting temperature from the 
fits of absolute yields is slightly higher than in case of ratios. This, as well as a poorer 
X^/^(i/=15.3/10, is mainly due to an apparent experimental upward bias of the yield of 
the S (and, to a lesser extent. A) hyperons. Without S's in the fit: T=160 MeV, yUfe=32 
MeV, y=2200 fm^, xVArd/=4.3/8. 

The main results for the best fits of absolute hadron yields are summarized in Table 5. 
The upper part in Table 5 is for the fits using mid-rapidity data, the lower part is for data 
integrated over Att. The fits of yields give similar temperatures as for ratios, but larger 
values for baryo-chemical potentials and larger values. One notices that, in general, 
X^/Nfif is close to or below unity, being significantly above unity only for the SPS energies, 
as for the case using ratios of yields. For all energies, the values are larger than in case 
of the fits of ratios. This is expected, as part of the discrepancy between data and model 
is cancelled in the ratios, for instance in case of antiparticle/particle ratios. These results 
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Fig. 23. Hadron yields with best fit at ^/sJfJ^=130 GeV. The dashed lines are for the best fit 
excluding the H hyperons. The yield includes both and 0. 

demonstrate that the use of ratios is less prone to biases. This is one of the reasons why 
we prefer to fit particle ratios rather than yields. 



B. Inclusion of 75 into the fits 

In Table 6, we present results obtained with mid-rapidity and 47r-integrated yields with 
a thermal model including a strangeness suppression factor 75. The results for the mid- 
rapidity data demonstrate that the key quantities, namely T and ^b, are very little influ- 
enced by whether or not 75 is included as a fit parameter. Also, the quality of the fits is 
not significantly improved by including 75. 

The fits of An yields including 75 leads to a substantial improvement of the fit quality 
only for the higher SPS energies, but still the fit seldom attains a good quality. For these 
energies, the resulting T (and also fib) values are larger than in the case without 75 
(Table 5). This effect is negligible for the lowest SPS energy and increases gradually as 
a function of energy, see Table 6. This is illustrated in Fig. 24 for the lowest and highest 
SPS beam energies. Contrary to expectations, the volume and 75 show a non-monotonic 
behavior as a function of energy (also apparent in the analysis of Becattini et al. [20]). In 
view of this, we regard V and 75 as compensatory parameters for an inadequate fit of Air 
data. As a consequence, the fact that the temperatures from the fits of in data using 75 
are very similar to those extracted from fits of mid-rapidity data should be viewed as a 
coincidence. 



Table 5 



Summary of the results of the thermal fits of yields. The upper part shows the results obtained 
using mid-rapidity data, the lower part is for data integrated over Air. 
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C. Importance of feed-down from weakly decaying states 



In Fig. 25 we show, as a function of energy, the contribution of weak decays to the yields 
of pions, protons and A hyperons calculated with the parametrizations for T and fib (Eqs. 
6 and 7 in Section 4) . Shown is the fraction of the total yield of these particles originating 
from weak decays. This fraction reaches the asymptotic values of 15%, 25% and 35% for 
TT, A and p, respectively, being significantly larger in case of antiparticles at AGS and SPS 
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Fig. 24. Hadron yields integrated over Air with best fits with (continuous hnes) and without 
(dashed hnes) 7s at the SPS energies of 20 and 158 AGeV. 

energies. Since this contribution is sizeable, a consistent comparison between data and 
model is required. 
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